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ABSTRACT 

Heat  transfer  rates  for  laminar,  convective  heat  transfer  in  the 
entrance  region  between  parallel  plates  are  investigated „    The  hydro- 
dynamic  solution  due  to  Bodia  /_2_/  was  used  in  the  solution  of  the  energy 
equation  in  finite  difference  form  on  a  digital  computer.    The  thermal 
boundary  conditions  include:    constant  heat  input,  constant  wall  tempera- 
ture, one  wall  constant  temperature  and  one  wall  insulated,  and  constant 
but  different  wall  temperatures  on  the  upper  and  lower  walls. 

The  approximate,  integral  methods  of  Siegel  and  Sparrow  /_\J ,  /_!_/ 
produce  results  that  are  in  close  agreement  with  the  solutions  in  this 
analysis  for  the  constant  heat  input  and  constant  wall  temperature  cases. 

The  scope  of  the  finite  difference  solution  is  limited  to  a  narrow 
range  of  Prandtl  numbers  near  unity,  due  to  the  small  grid  sizes  required 
for  convergence  at  small  Prandtl  numbers  and  to  the  overly  low  transfer 
rates  indicated  near  the  entrance  for  high  Prandtl  numbers,  which  is  a 
result  of  the  "finite  starting  length," 
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I.     INTRODUCTION 

Efficient  design  of  compact  heat  exchangers  such  as  are  employed 
for  gas  turbine  regenerators,  among  other  applications,  requires  maximum 
utilization  of  the  high  transfer  rates  available  in  the  entrance  region.    In 
this  region  the  velocity  and  temperature  profiles  of  the  entering  fluid  are 
undergoing  rapid  transition  from  their  uniform  distribution  at  entrance  to  the 
fully  established  distributions  encountered  farther  downstream.    The  well 
known  Graetz  solution,  in  which  only  the  temperature  development  is  con- 
sidered, does  not  yield  adequate  results  when  the  velocity  and  temperature 
distributions  are  simultaneously  developing.    This  is  particularly  true  if  the 
entrance  region  represents  a  considerable  percentage  of  overall  length  and 
the  two  profiles  develop  at  approximately  the  same  rates.    This  situation  is 
frequently  encountered  with  gases  where  the  Prandtl  number  is  approximately 
one. 

The  purpose  of  the  present  analysis  is  to  present  a  finite  difference 
solution  for  the  heat  transfer  rates  between  parallel  flat  plates  for  several 
different  boundary  conditions  corresponding  to  various  constant  wall 
temperatures  and  constant  heat  input  configurations.    The  approach  to  the 
problem  was  suggested  by  Miller  /_!_/ ' ,  which  is  adapted  from  the  numerical 
procedures  applied  to  the  entrance  region  of  the  circular  tubes  by  Kays  /_2J \ 


While  an  analysis  such  as  the  present  one  can  hardly  compete  for 
simplicity  with  approximate  treatments  such  as  the  Graetz  analysis ,  or 
the  several  energy  integral  analyses  commonly  employed  in  engineering 
design,  it  does  provide  an  exact  solution  for  a  limited  number  of  cases 
which  serves  as  a  standard  with  which  the  accuracy  of  the  various  approxi- 
mate methods  may  be  assessed. 


II.    ANALYSIS 

A.         Governing  Equations 

The  governing  equations  are  the  energy,  momentum  (Navier-Stokes) 
and  continuity  equations .    In  order  to  reduce  the  complexity  and  coupling 
of  the  equations,  the  following  assumptions  are  made. 
The  flow  is  assumed  to  be: 

1„    steady,  two-dimensional 

2 .  laminar 

3.  incompressible 

In  addition  it  is  assumed  that: 

1.  thermal  diffusivity  (o<-  )  is  constant 

2.  convective  heat  transfer  is  large  compared  to  radiation,  axial 
conduction,  and  viscous  dissipation. 

The  governing  equations  in  reduced  form  may  then  be  expressed  as: 
Energy: 

U  *£.     +  -V  il      r    CC      «  (1) 

Momentum: 

ul%         Ti         e     a*  *«f  (2) 


Continuity: 


iii  +    ££    =  0 


(3) 


Due  to  the  assumption  of  incompressibility  the  above  equations  are 
no  longer  coupled.    The  temperature  development  remains  dependent  upon 
the  velocity  profile,  although  dependence  of  velocity  on  temperature  has 
been  removed . 

The  boundary  conditions  associated  with  equations  (1),   (2),  and 
(3)  for  flow  between  parallel  flat  plates  with  uniform  velocity  and  tempera- 
ture at  entrance  are: 

1 .  At  the  entrance  (  x  =  0) 

t  =  t    =  C 
°  1 

u  =  uQ  =  C2 
v  =  vQ  =  0 

2 .  At  the  walls  (  y  =  +  d) 


v  =  vw  =  0 


u  =  uw  =  0 


t  =  tw 


3.    On  centerline 


v  =  vc  =  0 


4  .    Applied  thermal  boundary  conditions 

(l)  symmetric  cases  (with  respect  to  duct  centerline) 
a.    constant  wall  temperature 


b.    constant  heat  input 
f-*H  =C4 

(2)  asymmetric  cases 
a.    one  wall  constant  temperature,  one  wall  insulated 

X    =    C5 


bo    constant  wall  temperature 

=    C 


f&Y 


6 


*w2     ~    °7 

where:    C_^  C7 
6         / 

Solution  of  the  energy  equation  requires  that  a  hydrodynamic 
solution  be  first  obtained . 

B.         Hydrodynamic  Solution 

Several  hydrodynamic  solutions  for  the  entrance  region  are  available. 
Solutions  have  been  presented  by  Bodoia  /_Zj ,  Schlichting  /  A_/ ,  Schiller  /5_/, 
and  Han  /6_/.    The  approximate  solutions  of  Schiller  and  Han  have  the 
advantage  of  being  analytic,  and  therefore,  less  cumbersome  to  use.    They 
are,  however,  inherently  less  accurate  than  the  numerical  solutions  presented 
by  Bodoia  and  Schlichting. 


Bodoia's  solution  was  obtained  by  evaluation  of  the  Prandtl  boundary 
layer  equations  in  finite-difference  form.    Velocity  profiles  were  obtained 
at  given  streamwise  locations  by  utilizing  matrix  methods  to  simultaneously 
solve  a  column  array  of  momentum  expressions  normal  to  the  flow. 

A  comparison  of  Schlichting's  series  solution  and  the  Bodoia  solution 
indicates  that,  in  addition  to  the  velocity  gradient  discontinuity  in  the 
Schlichting  solution  where  the  upstream  and  downstream  solutions  are  joined, 
Schlichting's  representation  of  the  pressure  gradient,  based  only  on  the 
centerline  velocity,  results  in  a  core  velocity  gradient  which  is  too  large 
and  a  wall  velocity  gradient  which  is  too  small  /3/. 

Comparisons  of  the  dimensionless  pressure  decrement  versus  length 
due  to  Bodoia,  Schlichting  and  Han  are  presented  in  Fig.  2.    Velocity  dis- 
tribution comparisons  between  Bodoia  and  Schlichting,  and  Schlichting 
and  Schiller,  are  shown  in  Figs  .  3  and  4. 

The  Bodoia  solution  was  adopted  for  the  present  analysis. 

C .        Energy  Solution 

Equations  (l)  and  (3)  may  be  written  in  dimensionless  form  by 

introducing  the  following  dimensionless  variables: 

T  =  t/t  X  =    ^,  x =    x/d  u  =  u/un 

Y  =  y/d  V  =  v/uQ 


In  dimensionless  form,  the  equations  are  then: 
Energy: 

v  *T+Rev£L  -.  J_  S,  ... 


iJURei}£=0  (5) 


Continuity: 
\ 

The  following  finite  difference  approximations  are  then  introduced: 


il-     AT     ..      Tx*i,y  ~  Tv 
£X        AX  AX 


frT  _     AT     _      Tx.y+i  "  Txky-i 
*y  "    AY  2  ay 


(6) 


(7) 


*X*  AY2  (8) 


$5Lz     &V    z      I/«*i,Y  -U«-i,v 
^X         AX  2  AX 


(9) 


iY=  4v    s     V/.y-Vx.Y-i  (10) 

^Y        AY  AY 


Note  that  (6)  is  a  forward  difference  and  (10)  is  a  backward 
difference,  whereas  the  remaining  equations  are  of  the  central  difference 


type.    These  have  been  expressed  in  this  manher  for  simplicity  and 
convenience.    The  larger  truncation  error  normally  associated  with  the 
forward  or  backward  difference  is  not  significant  here  since  practical  grid 
dimensions  require  that: 

AY   »  AX 

and ,  due  to  the  nature  of  the  flow: 

AU    \\     AV 
AX  AY 

Introducing  (6)  through  (10)  into  (4)  and  (5)  and  solving  for  the 

appropriate  term  yields: 

T  -  f  AX  AX  ReV]    T  , 

'XU.Y    -  L   UPRAY2  VZLYJ        X'Y'X 


Vv  -Vw+  4r[u^-H 


(12) 


which  constitute  the  computing  equations . 

Note  that  the  Reynolds  number  (R  ) ,  although  appearing  in  the 
above  equations  „  remains  only  for  convenience  and  could  be  eliminated 
It  is  no  longer  a  parameter  in  the  set  of  equations  having  been  included 
in  the  dimensionless  length  variable  (X)  . 


A  mixed  mean  temperature  (t    )  may  be  expressed  in  the  following 

way: 

/Co  p  U t  dA 

tm=      V        (13) 

m         /cp^Uc/A 

Removing  the  constants  from  within  the  integral,  and  simplifying,  a 

dimensionless  mixed  mean  temperature  may  be  written  in  finite  difference 

form: 

Y--1 

Tm  =    -J£2 (14) 

Z    Ux,Y 
Y=o 


Due  to  the  no-slip  condition  imposed  upon  viscous  flows,  heat 
transfer  to  or  from  the  fluid  at  the  wall  may  be  expressed  in  terms  of 
conduction  through  a  thin  film  of  stagnant  fluid.    Thus: 

r 

q'-«lA-k(^)Jrt  (15) 

which  must  equal  the  heat  flux  due  to  convection, 


q     =  q/A  -  h(tw  -  tm)  (16) 

Defining  the  local  Nusselt  number  as; 

NUx=     ^  (17) 

Combining  (14)  and  (15) , 

(a) 


x 


t»-t 


w      fcm 


A  Nusselt  number  based  on  the  hydraulic  diameter  (d)  instead  of 
half  the  gap  width  (d)  in  finite  difference  form  becomes; 

N        =       — ^~ (IS) 

UX  '  w  \rr\ 


The  hydraulic  diameter  was  introduced  in  (19)  to  enable  comparisons 
with  previous  solutions.    The  conversion  is: 
D  =  4d 

In  addition,  for  the  asymmetric  constant  wall  temperature  case,  it 
is  convenient  to  use  a  dimensionless  parameter  other  than  the  Nusselt 
number.    A  dimensionless  heat  flux  parameter  may  be  defined  by  using 
(15)  above: 


0 


kt 


d     Uy/Y-o 


=  4  ("5yJy«o 


D.        Thermal  Boundary  Conditions  In  Finite  Difference  Form 

Constant  Wall  Temperature femmetric  Case)    Since  the  temperature 
distribution  is  symmetrical  with  respect  to  the  centerline  ,  computations 
were  performed  in  only  the  upper  half  of  the  duct.    A  coordinate  system  was 
adapted  with  the  X  axis  on  the  centerline  of  the  channel.    The  thermal  boundary 
condition  at  the  wall  was  applied  as: 


10 


t    /tn    =    T„    =    N 
w    o  w 

where  N  is  a  positive  integer. 

Constant  Heat  Input  (jfyjnmftttfis Case)    With  the  axis  system  as  above , 

the  slppe  of  the  thermal  profile  was  maintained  constant  at  the  wall  by 

adding  a  constant  to  the  calculated  temperature  one  step  from  the  wall, 

thus: 

T     =  T     -        +  M 
lw      lw     Ay 

where:      M    =    (  AT/AYV   AX 

W 

Constant  Wall  Temperature  with  One  Wall  Insulated .    For  this  case 
the  coordinate  system  was  redefined  so  thfe  X  axis  coincided  with  the  lower 
wall.    The  lower  wall  boundary  condition  was  taken  to  be: 

t      /t    =  T        =  N 

wy    o        wj_ 

while  zero  slope  at  the  upper  wall  was  obtained  by  equating  the  wall 
temperature  to  the  calculated  temperature  one  step  from  the  wall: 

tw2//t0  =  Tw2  =  T&2  -  AY  : 

Constant  Wall  Temperaturei(Agymmfftrjb  gajse),.  With  the  axis  system 
the  same  as  for  the  insulated  wall  case,  thermal  boundary  conditions  were 
established  as; 

tw/to  =  Tw1    =N1 

VA  =  TW2    =N2 

where  N-,  and  N2  are  positive  integers  and  N,    4  N2- 

11 


Noting  that: 

twi  _  t0  =        twiAo  "  1      =      Nx  -  1 


tW2  -  tQ  tW2/tQ  _  i  N2  -  1 

Values  of  N1  -  1/  N2  -  1  of  1/2 ,  1/3,  1/4  were  investigated. 

E.         Computational  Procedure 

The  finite  difference  solution  was  obtained  by  using  a  CDC  1604 
digital  computer  to  solve  the  equations  at  intervals  of  AX  =  CL001  and 
AY  =  0.1 .    Velocity  distributions  were  obtained  by  interpolating  plotted 
curves  of  Bodoia's  tabulated  results. 

Since  the  entire  solution  is  dependent  upon  the  first  calculations 
at  the  entrance,  starting  values  were  obtained  by  reducing  the  grid  size 
to  AX  =  0.0001  and  A  Y  =  0.05  for  the  first  20  streamwise  stations  from 
X  =  0  to  0,0020. 

The  continuity  equation  (2)  was  evaluated    and  the  cross-wise 
velocity  component  (V)  was  introduced  into  the  energy  equation  (11)  at 
each  grid  point.    The  energy  solution  was  then  obtained  at  intervals  of 
Y  between  centerline  and  the  wall  at  each  streamwise  station. 

Solutions  for  each  thermal  boundary  condition  were  obtained  for 
Prandtl  numbers  of:    0.5,  0.7,  1.0,  1.6,  3.2  and  10. 

Solutions  for  Prandtl  numbers  less  than  0.5  are  excluded  by  the  grid 
size.    This  may  be  seen  from  inspection  of  the  coefficient  of  the  second 
temperature  term  in  (ll)#  which  indicates  the  interdependence  of  PR, AX, 

12 


and  AY.    In  order  to  obtain  convergence,  PR  must  be  large  enough  to 
prevent  this  term  from  becoming  negative.    Thus  low  Prandtl  numbers, 
approaching  those  of  liquid  metals,  cannot  be  effectively  handled  by 
finite-difference  methods  owing  to  the  large  number  of  calculations 
required  by  the  small  grid. 

The  computer  program,  written  in  Fortran  60,  is  included  in  the 
Appendix  „ 
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III.    RESULTS 

Tables  I  through  VI  contain  computed  and  extrapolated  values  of 
the  heat  transfer  rates  as  a  function  of  distance  from  the  entrance  for  each 
of  the  thermal  boundary  conditions  for  a  Prandtl  number  of  0.7  „    Figures 
5  through  29  contain  curves  of  the  transfer  rate  versus  downstream  position 
and  temperature  and  velocity  profiles  for  Prandtl  numbers  of  0.5,  0,7,  1.0, 
1.6,  3o2  and  10. 

With  the  exception  of  the  Asymmetric  Constant  Wall  Temperature 
Case,  which  is  expressed  in  terms  of  the  dimensionless  heat  flux  parameter, 

* 

q  ;  the  heat  transfer  rates  are  presented  as  local  Nusselt  number  variation 
with  the  length  parameter,  L.    Whereas,  the  computations  were  carried  out 
using  the  dimensionless  length  parameter  X,  which  is  normally  employed 
for  hydrodynamic  developments;  L,  based  on  the  hydraulic  diameter  is 
customarily  employed  in  heat  transfer.    The  conversion  is: 
X  =  16L 
The  use  of  all  finite  starting  length"  in  the  finite  difference  equations 
results  in  a  finite  initial  value  for  the  Nusselt  number  at  X  =  0  rather  than 
the  theoretically  infinite  value  at  the  entrance.    This  causes  a  nearly  flat 
slope  near  the  entrance  and  an  artificial  inflection  point  in  the  Nusselt 
number  versus  length  curve.    The  inflection  point  occurs  progressively 
farther  downstream  with  the  increasing  Prandtl  number,  thereby  limiting  the 
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the  useful  span  of  the  curves.    A  reduction  of  grid  size  from  AX  =  .001 , 
AY  =  .  1  to  4  X  =  . 000 1 , 4 Y  =  .05  was  observed  to  double  the  initial 
Nusselt  number  from  40  to  80,  arid  move  the  inflection  point  upstream. 
This  increased  the  useful  span.    However,  the  approach  obviously  reaches 
a  practical  limit  very  rapidly  due  to  the  number  of  calculations  and  the 
accuracy  required  of  the  hydrodynamic  solution  very  close  to  the  entrance. 

In  general,  for  the  three  thermal  boundary  conditions  in  which  the 
Nusselt  number  was  used  as  a  heat  transfer  parameter,  the  curves  for 
Prandtl  numbers  of  0.5,  0.7,  lo0  and  1,6  behave  as  expected  for  values 
of  L>  10"    .    For  Prandtl  numbers  of  3.2  and  10,  L  =  10"3  is  the  approxi- 
mate lower  limit  of  validity. 

Constant  Heat  Input .    A  comparison  of  local  Nusselt  numbers  obtained 
in  the  present  analysis  with  those  reported  by  Siegel  and  Sparrow  /_!_/,  and 
Han  /6/,  contained  in  Fig.  5,  indicate  close  agreement  with  the  results  of 
Siegel  and  Sparrow.    Their  solution  employed  a  simplified  energy  equation 
in  which  the  temperature  distribution  in  the  boundary  layer  was  expressed 
as  a  series  of  polynomials  in  the  transverse  coordinate  using  the  down- 
stream station  as  a  parameter.    The  hydrodynamic  Solution  employed 
Schiller's  approximation. 

The  comparison  with  Han's  integro -numerical  solution  is  not  as  close 
as  might  appear  in  Fig.  5  since  this  curve  is  for  a  Prandtl  number  of  0.8 
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and  should  lie  above  the  other  two  curves  which  are  for  a  Prandtl  number 
of  0.7. 

The  computed  curves  of  local  Nusselt  number  versus  length  contained 
in  Fig .  6 .  for  Prandtl  numbers  of  0.5,  0.7,  1.0,  1.6,  3  o  2 ,  and  1 0 , i may  be 
approximated  within  10  per  cent ,  for  0.5^  PR  ^  1.6,  and  L  >  1 0~     by  the 
following; 


N„      =    8.24    +  -0186 


(PR/L) 


Ux  1+  .0178  (PR/I.)'6 

Velocity  and  temperature  profiles  for  constant  heat  input  are  contained 
in  Figs.  7  through  10. 

Constant  Wall  Temperature  (Symmetric  Case) .    Local  values  of  Nusselt 
numbers  for  this  boundary  condition  are  compared  in  Fig .  1 1  with  those  of 
Sparrow/ 5/,  which  were  computed  using  the  Karman-Pohlausen  method  and 
Schiller's  hydrodynamic  solution.    Again  close  correlation  is  indicated 
between  the  two  results. 

The  local  Nusselt  number  curves  of  Fig.  12  may  be  approximately 

—3 

represented  within  JO  per  cent  for    0 .  5  £  PR  i  1 . 6  ,  and  L  >  1 0  '    by  the 

following; 

.0148  (PRA) 


N„      =    7„60    + 


ux  1  +  .0143  (PR/L)-6 

Velocity  and  temperature  profiles  for  this  boundary  condition  are 
contained  in  Figs.  13  through  16. 
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One  Wall  Constant  Temperature  -  One  Wall  Insulated,    The  results 

of  this  computation  are  presented  in  Fig,  17.    The  local  Nusselt  number  may 

—S 
be  approximated  within  10  per  cent  for  0.5^  PR-  1.6,  and    L  >  10       by; 

N        =    4   84    +        .0155  (  PR//L) 
U*  1    +     .012  (PR/L)'b 

The  velocity  and  temperature  distributions  are  illustrated  in  Figs. 
18  through  21  „ 

Constant  Wall  Temperature  (Asymmetric  Case).    Dimensionless  heat 

flux  values  for  this  boundary  condition  are  presented  in  Figs.  22,  23  and  24 

for    Ni  -  1       ratios  of  1/2,  1/3,  and  1/4. 
H2-l 

Note  that  the  upper  (hotter)  wall  curve  is  smoothly  asymptotic  to 
the  value  predicted  for  the  steady  state  condition.    The  lower  curve 
exhibits  a  much  shallower  initial  slope,  passes  through  an  inflection 
point,  and  then,  approaches  a  negative  asymptotic  value  of  the  same 
magnitude  and  at  the  same  rate  as  the  upper  wave. 

Velocity  and  temperature  profiles  for    N\     •  1  =    1/4  are  presented 
in  Figs.  25  through  28. 


N2-l 


The  behavior  of  the  local  Nusselt  number  for  this  asymmetrical  case 
is  depicted  in  Fig.  29.    The  lower  curve,  representing  the  local  Nusselt 
number  variation  for  the  lower  (cooler)  wall,  does  not  depict  the  actual 
physical  situation. 
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Both  the  denominator  and  the  numerator  of; 


4(  Tw-Tw-^y) 

"%"  T*-T.  (18) 


approach  zero  at  different  rates.    The  curve  tends  to  positive  infinity  as 
the  denominator  approaches  zero,  returns  from  negative  infinity,  and 
becomes  zero  as  the  numerator  goes  to  zero;  and  finally  approaches  the 
actual  steady  asymptotic  value. 
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IV.     CONCLUSIONS 

From  the  results  of  this  study  it  may  be  concluded  that: 
1.      The  finite  difference  method,  as  employed  in  this  analysis, 
yields  results  which  are  consistent  with  previous  solutions 
for  other  geometries  /_2j  and  which  are  asymptotic  to  the 
known  fully  developed  values.    The  method  is,  however, 
limited  to  a  small  range  of  Prandtl  numbers  near  unity  by  the 
small  grid  sizes  necessary  outside  this  range »    For  small 
Prandtl  numbers  the  grid  must  be  very  small  to  obtain  con- 
vergence,  while  for  large  Prandtl  numbers  practical  grid 
dimensions  result  in  significant  inaccuracies  near  the 
entrance  due  to  the  "finite  starting  length". 
2  .      The  results  of  the  approximate  integral  methods  suggested 
by  Siegel  and  Sparrow  /_!_/  for  constant  heat  input,  and  by 
Sparrow  /5_/  for  constant  wall  temperature,  compare  very 
favorably  with  the  finite  difference  solutions  reported  here. 
This  additional  substantiation  of  the  methods  suggested  by 
these  authors  is  significant  because  of  the  practical  utility  of 
their  approach,which  requires  relatively  few  calculations,  is 
capable  of  handling  a  wide  range  of  Prandtl  numbers  ,  and  can 
predict  heat  transfer  rates  very  close  to  the  entrance. 
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TABLE  I 
CONSTANT  HEAT  INPUT 


PR  =  0 . 7 

X 

L 

Nu 
ux 

T 
w 

0016 

. ,0001 

46„2 

lo09 

0032 

.0002 

30.0 

1.14 

0064 

.0004 

23.7 

1.18 

0128 

.0008 

18.0 

1.24 

0160 

.0010 

16.6 

K26 

0320 

.0020 

12.9 

K35 

0640 

.0040 

10„3 

1.47 

1280 

.0080 

8.87 

1.62 

1600 

.0100 

8.62 

1„67 

Asymptotic  to     8.24 
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TABLE  II 
CONSTANT  WALL  TEMPERATURE 
(SYMMETRIC  CASE) 
PR  =  0,7 

X  L  Nv 

0016  .0001  34. 


ux 


.0032  .0002  26„0 

.0064  .0004  1909 

.0128  ,0008  15.3 

.0160  .0010  13.7 

.0320  .0020  10.5 

.0640  .0040  8.58 

.1280  .0080  7.78 

.1600  .0010  7.63 

Asymptotic  to   7.60 
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TABLE  III 
ONE  WALL  CONSTANT  TEMPERATURE  -  ONE  WALL  INSULATED 


PR  =  0 . 7 

X 

L 

ux 

T 
w 

.0016 

.0001 

33.0 

1.000 

.0032 

.0002 

24,0 

1.000 

.0064 

.0004 

19.2 

1.000 

.0128 

.0008 

13.0 

1.000 

.0160 

.0010 

11.6 

1.000 

.0320 

.0020 

9.05 

1.000 

.0640 

.0040 

7.20 

1.000 

.1280 

.0080 

5.93 

1.001 

.1600 

,0100 

5.55 

1.021 

.4000 

.0250 

4.96 

1.139 

Asympl 

:otlc  to-    4.84 

2.000 
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X 
.0032 
00064 
.0128 
.0160 
.0320 
.0460 
.1280 
.1600 
o3200 
.6400 
.9600 


TABLE  IV 
CONSTANT  WALL  TEMPERATURE 
(ASYMMETRIC  CASE) 


*U  -  1  =  1/ 

2    PR  =  0  . 7 

N2  -  1 

L 

* 
'1 

* 

q2 

.0002 

26.2 

57.5 

.0004 

17.3 

41.4 

.0008 

11.9 

28.7 

.0010 

10.7 

25.9 

.0020 

8.08 

18.7 

.0040 

6.18 

13,7 

.0080 

4.73 

10.1 

.0100 

4.24 

9.15 

.0200 

2.24 

6,35 

.0400 

-.19 

3,81 

.0600 

-1.23 

2,76 

Asymptol 

tic  to    -2.00 

+  2.00 
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TABLE  V 

CONSTANT  WALL  TEMPERATURE 

(ASYMMETRIC  CASE) 

Nl    -     1      =1/3  PR  =  0.7 
N2    -    1 

X                                              L  qx* 

.0032                                      .0002  28.0 

.CM364                                       .0004  17.3 

.0128                                     .0008  11.8 

.0160                                     .0010  10.8 

.0320                                     .0020  8.09 

.0640                                     .0040  6.18 

.L280                                     .0080  4.69 

.1600                                     .0100  4.15 

.3200                                     .0200  1.65 

.6400                                     .0400  -1.58 

.9600                                     .0600  -2.97 

Asymptotic  to  -4.00 


* 
q2 

86 

.0 

62 

.0 

43 

.0 

38 

.9 

28 

.1 

20 

.6 

15 

2 

13, 

8 

9, 

84 

6 

41 

5. 

02 

+4 

00 

25 


TABLE  VI 

CONSTANT  WALL  TEMPERATURE 

(ASYMMETRIC  CASE) 


N 


1    -    1  =  1/4  PR    =   0 . 7 


N2    -    1 

X  L 

,0032  .0002 

,0064  .0004 

0128  .0008 

0160  .0010 

,0320  .0020 

0640  .0040 

1280  .0080 

1600  .0100 

3200  .0200 

6400  .0400 

9600  .0600 


Asymptotic  to 


* 

ql 

* 
q2 

28.2 

116.0 

16.8 

82.9 

12.2 

57,7 

10.5 

50.4 

8.09 

37.5 

6.17 

27.5 

4.66 

20.3 

4.06 

18.4 

1.06 

13.3 

-2.96 

9.03 

-4.71 

7.27 

-6.00 

+6.00 
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TEMPERATURE     AND  VELOCITY     PROFILES 

CONSTANT    HEAT    INPUT 
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TEMPERATURE     AND    VELOCITY    PROFILES 
X=  0.050  L  =  0.003125 
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TEMPERATURE    AND    VELOCITY      PROFILES 
,       X=  0.250  L  =0.0156 
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FIGURE     10 

TEMPERATURE  AND  VELOCITY     PROFILES 

X=  1.00  L  =0.0625 
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TEMPERATURE    AND   VELOCITY      PROFILES 
CONSTANT     WALL     TEMPERATURE 

X=0.005  L  =0.0003125 
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TEMPERATURE    AND   VELOCITY      PROFILES 
CONSTANT    WALL    TEMPERATURE 
X=  0.250  L  =0.0156 
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TEMPERATURE     AND    VELOCITY     PROFILES 

CONSTANT     WALL    TEMPERATURE 
i     X=  LOO  L  =0.0625 
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FIGURE    18 

TEMPERATURE    AND   VELOCITY    PROFILES 
ONE   WALL    CONSTANT   TEMPERATURE 
ONE    WALL    INSULATED 
X=  0.005  L  =  0  .0003125 
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TEMPERATURE    AND    VELOCITY     PROFILES 
ONE   WALL  CONSTANT    TEMPERATURE 
ONE    WALL    INSULATED 
X=0.050  L  =0.003125 
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TEMPERATURE     AND    VELOCITY      PROFILES 
ONE    WALL    CONSTANT    TEMPERATURE 

ONE     WALL    INSULATED 
X- 0.250  L  =0.0156 
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FIGURE      21 

TEMPERATURE     AND  VELOCITY     PROFILES 
ONE    WALL    CONSTANT    TEMPERATURE 
ONE    WALL   INSULATED 
X=  1.00  L  =0.0625 
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FIGURE      25 

TEMPERATURE   AND  VELOCITY     PROFILES 
CONSTANT    WALL   TEMPERATURE 
(ASYMMETRIC    CASE) 

X=0.005  L  =0.0003125 
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FIGURE    26 

TEMPERATURE     AND   VELOCITY      PROFILES 

CONSTANT    WALL    TEMPERATURE 

(ASYMMETRIC   CASE) 

X=  0.050  L  =0.003125 
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FIGURE    27 

TEMPERATURE    AND    VELOCITY     PROFILES 
CONSTANT    WALL    TEMPERATURE 

(ASYMMETRIC  CASE) 

X=  0.250  L    =0.0156 
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FIGURE     28 

TEMPERATURE    AND   VELOCITY      PROFILES 
CONSTANT     WALL     TEMPERATURE 
(ASYMMETRIC    CASE) 

X=  1.00  L=  0.0625 
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APPENDIX  A 
SAMPLE  FORTRAN  60  PROGRAM 
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/ 


..JOB0191F,LUNDBERG 

PROGRAM  HEATX23 
C      CONSTANT  HEAT  INPUT 
C      BODOIA  VELOCITY  DATA 

C      THIS  PROGRAM  USES   A  REDUCED  GRID   OF.  X  =  .000 1 , Y  =  . 05 

C      FOR  THE  FIRST  20  STREAMWISE  CALCULATIONS. 

C      GRAPH  OUTPUTS  OF  NUSSELT  NUMBER  VERSUS  LENGTH  PARAMETER 

C      FOR'PRANDTL  NUMBERS  OF  0.  5  ,0.7 » 1.0 ,1 .6 , 3 .2 ,AND  10»   '■ 

C      ARE  OBTAINED  TO  50  PER  CENT  OF  THE  HYDRODYNAMIC 

C      DEVELOPMENT.   PRINT  OUTPUT  INCLUDES  NUSSELT  NUMBER, 

C      SOPE  OF  THE  TEMPERATURE  PROFILE  AT  THE  WALL  ,  AND  THE 

C  .    SLOPE  OF  THE  TEMPERATURE  PROFILE  AT  THE  WALL, AND  THE  MEAN 

C      MIXED  MEAN  TEMPERATURE.   IN  ADDITION,  VELOCITY  AND 

C      TEMPERATURE  VALUES  ARE  PRINTED  AT  EACH  GRID  POINT  AT 

.C      SELECTED  STREAMWISE  LOCATIONS. 

DIMENSION  TOO, 23)  ,U(30,23)  ,V(30,23)  ,Y(23)  ,X(500) 
1,XN(500) ,ITITLE(12) 

READ  500, ( ITITLE( I )»  1=1,6) 

READ  500, ( ITIJLE( I ),  1=7,12) 
500  FORMAT  (6A8) 

DO  750  KK=1,6 

GO  T0(40, 41,42,43, 44, 45), KK 

40  PR=0.7 
GO  TO  46 

41  PR=0.5 

GO  TO  888 

GO  TO  46  -  • 

42  PR=1.0 
GO  TO  46 

43  PR=1.6 
GO  TO  46 

44  PR=3.2 
GO  TO  46 

45  PR=10. 

46  PRINT  523 
523  FORMAT  (1H1) 

PRINT  17,  PR 
17  FORMAT  (/////10X,17H  PRANDTL  NUMBER  =,F5.3) 
L=l 
GO  TO  150 

700  READ  701,  ((U(I,J),  J=13»22)»  1=2,22) 

701  FORMAT   (10F6.4) 
DO  710  1=2,22 

DO  711  J=l,12 
711  U( I ,J)=U( 1,13)  ' 
710  CONTINUE 

DELX=.0001 

DELY=0.05 

R=1000. 

DO  110  L=2,21 
150  UTM=0.0 

UM=0.0  .      . 

XL  =  L 

X(L)=(XL-1.0)*0.0001 


60 


\ 


— -r- — irrrrrfin — rrrrvr 


i        j  t    l      i*i      i 


1 


108 
109 


113 


112 

114 

111 


131 
118 


119 


DO  111  J=3» 

XJ  =  J 

Y(.L)  =  (XJ-3. 

Y(23)=1.0 

U( 1,J)=1. 

U(l,23)=1.0 

V(1,J-1)=0. 

V( 1»JJ=0. 

T(1,J)=1.0 

TU,J+1)  =  1. 

T( 1,J-1)=1. 

T(L,1)=T(L» 

U(L»1)=U(L» 

U(L,23)=0.0 

U(L,2)=U(L, 

V(L»3)=0.0 

V(L,2)=V(L, 

THERMAL  BOU 

T(L»23)=T(L 

T(L»2)=T(L, 

IF(L-l)  109 

V(L,J)=  V(L 

A=  DELX/  (U 

V(L»23)=0.0 

B=  (DELX*R* 

.  T(L+1,J)=(A 
IFU-3)  113 
UTM=UTM+.5* 
UM=UM+.5*U( 
GO  TO  114 
UTM=UTM+U(L 
UM=UM+U(L,J 
TM=UTM/UM 
CONTINUE 
Tl=    (  T(L 
T2=     T<L, 
XN(L)=4.0*T 
PRINT  118, 
FORMAT  (// 

1F8.5,5X,8H 
PRINT119,  ( 
PRINT119t  ( 
FORMAT  (/9H 


22 

0)*0.05 


0 

0 

5)  ' 

5>  * 

4) 

4) 

NDARY  CONDITIONS 

,2.2)+. 05 

4) 

,109,108 

»J-4)  +  DELY/(2.*DELX*R)  *  ( U ( L-l , J ) -U( L+l . J ) ) 

(L,J)*PR*DELY**2) 

V(L,J))  /  (2.*DELY*U(L,J) ) 

-B)*T(L,J+1  )+(l.-2.*A)*T(L»J)+(A+B)*T(L,J-U 

♦113,112 

U(L,J)*T(L,J) 

L,J) 

,J)*T(L»J) 

)  ' 


,23)-T(L»22) )  /  DELY 

23)  -  TM 

1/T2 

X(L) ,XN(L) ,T2,T1,TM 


4H  X  =,F6.4,5X,14H 
SLOPE  =,F8.5,5X,5H 


NUSSELT  NO. 
TM  =,F8.5) 


=,F8.5,5X,8H  TW-TM  *, 


PRINT120, • ( 
PRINT120,  ( 

120  FORMAT  (/  9 
PRINT121,  ( 
PRINT121,  ( 

121  FORMAT  (/  9 
IF(L-l)  700 

110  CONTINUE 

CALL  DRAW  < 
DO  720  N*2» 
M=2*N-1 


T(L,J) 
T(L,J) 
TEMP. 
U(L,J) 
U(L,J) 
H  X  VEL 
V(L,J) 
V(L,J) 
H  Y  VELv 
,700,110 


J=l,12) 

J=13,23) 

=,12F7.4) 

J=l,12) 

J=13,23) 

=,12F7.4) 
J=l,12) 
J=13,23) 

=,12F7.4) 


21»X,XN,MOD,0,LAB,ITITLE,.04,10.,0,0,0,0,6,  6,1, LAST) 
11 


61 


i — r 


-■rr 


771 

720 

85 

70 


U(1»N)=U(11»M) 

U(2,N)=U(21,M) 
T(2,N)=T(21,M) 
PRINT  771,  Mi  U( 
FORMAT  (I10,2F10 
CONTINUE 
L  =  2 

DO  10  M=2t250 
XM  =  M 
Y=0.0 
DELY=0.1 
DELX=0.001 
X(M)=DELX*XM 
I=L  +  1 
N3N=M 
N  =  L  +  1 

IF(M-  99)  72,73, 

READ  1  (U( I ,11)  , 

1 ,U( I  ,5)  ,U( I ,4)  ,U 

1  FORMAT  (10F6.4) 

73  IF(M-99  )  50,51, 

IF(M-149)  51,53, 

IF(M-199)  53,55, 

UU=XM-  99. 


2,N) ,T(2,N) 
.•5) 


82 
72 


52 
54 
51 


53 


U(N,2) 

U(N,3) 

U(N,4) 

U(N,5) 

U(N,6) 

U(N,7) 

U(N,8)  = 

U(N,9)  = 

U(N,10)= 

U(N,11)= 

GO  TO  50 

UU=XM-149. 


=1.4388+. 
=1.4292+. 
=1.3993+. 
=1.3454+. 
=1.2628-. 
=1.1467-. 
=  .9938-. 
=  .8022-. 

.5720-. 

.3042-. 


U(N 
U(N 
U(N 
U(N 
U(N 
U(N 
U(N 
U(N 
U(N 
U(N 
60 


,2) 
,3) 
,4) 
,5) 
,6) 
,7) 
,8)  = 
,9)  = 
,10)  = 
,11)  = 
TO  50 


=1.4758+. 
=1.4629+. 
=1.42  39+. 
=1.3573+. 
=1.2611-. 
=1.1336-. 
=  .9733-. 
=  .7796-. 

.5526-. 

.2926-. 


55  UU=XM-199. 


U(.N,2) 
U(N,3) 
U(N,4) 
U(N,5) 
U(N,6) 
U(N,7) 


=1.4903+. 
=1.4762+. 
=1.4336+. 
=1.3619+. 
=1.2604-. 
=1.1284-. 


73 

U( 1,10)  ,U( 1,9) 

( I  ,3) ,U( 1,2)  ) 

52 
54 
55 

000740*UU 
000674*UU 
000492*UU 
000238*UU 
000034*UU 
000262*UU 
000410*UU 
000452*UU 
000388*UU 
000230*UU 


000290*UU 
000266*UU 
000194*UU 
000092*UU 
000014*UU 
000104*UU 
000160*UU 
000176*UU 
000150*UU 
000092*UU 


0000121*UU 
0000110*UU 
0000080*UU 
0000026*UU 
0000005*UU 
0000042*UU 


,U( I ,8) ,U( I ,7) ,U( I .6) 


62 


~C" 


-.  . ,  - ... 


U(N»8)  =  .9653-.0000066*UU 
U(N»9)  =  .7708-.0000072*UU 
U(.N.10)=  .5451-.0000063*UU 
U(N,11)=  .2880-.0000475*UU 
50  UTM=0.0 
UM=0.0 

00  11  J=2,ll 
XJ  =  J 

Y< J)«DELY*XJ-.2  '   ' 

U(L»12)=0.0 
V(L,2)=0.0 

V(L»1)=V( 1,3)  ' 

YQ2)  =  1.0 

T(L,12).=  T(L,11)+0.1 
T(L,1)=T(L,3) 
U(L,1)=U(L,3) 
8  V(L,J)=  V(LtJ-l)  +  DELY/(2.*DELX*R)  *  { U ( L-l , J ) -U ( L  +  l ,  J ) ) 
>=  DELX/  (U(L,J)*PR*DELY**2) 
V(L,12)=0.0  , 

B=  (DELX#R*V(L,J) )  /  ( 2 .*DELY*U ( Lt J ) ) 
.   T(L+1,J)=(A-B)*T(L,J+1 )  +  Q.-2.*A)*T(L»J)  +  <A+B)*T(L,J-l) 
IFU-2)  13,13,12 

13  UTM=UTM+.5*U(L,J)*T(L»J) 
UM=UM+.5*U(L»J> 

GO  TO  14 
12.  UTM=UTM+U(L,J)*T(L,J)  .> 

UM=UM+U(L,J) 

14  TM=UTM/UM 
11  CONTINUE 

Tl=   (  T(L,12)-T(L,11)  )  /  DELY 
T2=    T(L,12)-TM 
XN(M)=4.0*T1/T2 
IF(N3N-50  )  31,31,30 

30  IF( (N3N/25)*25-N3N)  15,31,31 

31  PRINT  18,X(M) ,XN(M>,T2,T1,TM 

18  FORMAT  (//  4H  X  =  ♦ F6. 3  ,5X  ,  14H  NUSSELT  NO.  =>F8.5,5X»8H  TW-TM  =  » 
1F8.5,5X,8H  SLOPE  =,F8.5,5X,5H  TM  s»F8.5> 

PRINT  19,  (T(L,J) ,  J=l,12) 

19  FORMAT  (/9H  TEMP.   =,12F7.4) 
PRINT  20,  (U(L,J) ,  J=l,12) 

20  FORMAT  (/  9H  X  VEL.  =,12F7.4) 
PRINT  21,  (V(L,J),  J=l,12) 

21  FORMAT  (/  9H  Y  VEL.  =»12F7.4) 

15  DO  60  J=2,ll 
T(L,J)=T(L+1,J) 
U(L-1,J)=U(L,J) 

60  U(L,J)=U(L+1»J) 
10  CONTINUE        '■ 

IF(KK-6)  760,761,761 

761  MOD=3 

GO  TO  762 
760  MOD=2 

762  CALLDRAW(497»X,XN»MOD,0,LAB,ITITLE,.04,10.,0,0,0»0»6»  6t  It  LAST). 
750  CONTINUE 

63 
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